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ABSTRACT 

We show that the smoothed particle hydrodynamics (SPH) method, used with individual time-steps in the 
way described in the literature, cannot handle strong explosion problems correctly. In the individual time- 
step scheme, particles determine their time-steps essentially from a local Courant condition. Thus they cannot 
respond to a strong shock, if the pre-shock timescale is too long compared to the shock timescale. This problem 
is not severe in SPH simulations of galaxy formation with a temperature cutoff in the cooling function at 10 4 K, 
while it is very dangerous for simulations in which the multiphase nature of the interstellar medium under 10 4 K 
is taken into account. A solution for this problem is to introduce a time-step limiter which reduces the time-step 
of a particle if it is too long compared to the time-steps of its neighbor particles. Thus this kind of time-step 
constraint is essential for the correct treatment of explosions in high-resolution SPH simulations with individual 
time-steps. 

Subject headings: galaxies:evolution — galaxiesTSM — galaxies:star formation — methods:numerical 
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1. INTRODUCTION 



Hierarchical (individual) time-step method ( e.g., McMillan 
U986tlHernquist & Katdll989HMakinolll991al) is widely used 
in simulations of galaxy formation and star formation based 
on smoothed particle hydrodynam ics (SPH) method dLucvl 
ll977tlGingold & Monaghanl 1977b . This method allows parti- 
cles to have different time-steps, and can significantly reduce 
the total calculation cost when there is a large variation in the 
timescales of particles. Almost all implementations of the in- 
dividual time-steps meth od used for particle system s violates 
Newton's third law (see iFarr & Bertschingerll2007l for one of 
exception). As long as physical quantities are integrated with 
sufficient accuracy, this violation is not a severe problem. 

However, it is not always possible to maintain the accuracy. 
To our knowledge, all existing implementations of individual 
time-step for SPH rely on the determination of the time-step 
at the end of previous time-steps. Therefore, if something 
unforeseen occurs during the time-step of one particle, the 
particle might fail to catch that event, resulting in a large in- 
tegration error. A supernova (SN) explosion is an example of 
such an event. A SN generates a small amount of very hot 
gas (T ~ 10 8 K) in a large clump of cold gas (T ~ 10 K). 
The difference in the time-steps of the hot gas and surround- 
ing cold gas particles becomes quite large (typically the dif- 
ference reaches ~ 10 3 ). Thus, hot gas particles step forward 
~ 10 3 or more time-steps before neighboring cold gas parti- 
cles respond to the SN event. This means that the evolution 
of both hot and cold gas particles is completely wrong, since 
the surrounding cold gas particles do not react the explosion 
for a duration much longer than the timescale in which the 
blast wave would propagate the inter-particle distance. This 
problem is not severe for ordinary SPH simulations of galaxy 
formation because of the temperature cutoff in cooling func- 
tions at 10 4 K, while it becomes very serious for simulations 
involving the multiphase nature of the interstellar medium un- 



Electronic address: saitoh.takayuki@nao.ac.jp, saitoh.takayuki@cfca.jp 

1 Center for Computational Astrophysics, National Astronomical Obser- 
vatory of Japan, Mitaka, Tokyo 181-8588, Japan 

- Division of Theoretical Astronomy, National Astronomical Observatory 
of Japan, 2-21-1 Osawa, Mitaka-shi, Tokyo 181-8588; and School of Phys- 
ical Sciences, Graduate University of Advanced Study (SOKENDAI) 



der 10 4 K, because the mach number can be very high. 

This problem of sudden change occurs in any dynamical 
simulation with individual time-steps, as long as the time- 
steps are determined with the usual explicit method. In princi- 
ple, a fully implicit me thod in which the tim e-step itself is also 
determined implicitly ( Makino ^t al.l2006|) or a method which 
satisfies Newton's third law (IFarr & Bertschingerl l2007l) can 
solve this problem, but there are no implementations of such 
methods for SPH simulations yet 3 . 

In simulation of star clusters, the SN and its kick introduces 
a sudden change in the orbit (and mass) of the exploded star. 
Here, a rather simple prescription in which either all stars or 
at least nearby stars are synchronized to the time of explosion 
and restart the integration has been used. This prescription, at 
least the version which synchronizes all particles in the sys- 
tem, is impractical for A/-body/SPH simulations of galaxies, 
because the number of SN events and therefore the increase 
in the calculation cost is too great. 

We propose a simple limiter for hydrodynamical time-steps 
in the individual time-step method which mitigates this prob- 
lem. With this limiter the behavior of an explosion integrated 
by individual time-steps becomes essentially the same as that 
integrated by global time-steps. In §2, we describe this lim- 
iter, and in §3 we report the result of numerical experiments. 

2. THE TIME-STEP LIMITER 

We denote the time-step of an i-th particle as dtj and that of 
a neighbor particle, with index j, as dtj. The basic idea of our 
limiter is to enforce the following conditions: 

dti < f dtj and, (1) 
dtj<fdt t , (2) 
where / is an adjustable parameter. We found / = 4 to give 
good results, from the perspective of the total energy and lin- 
ear momenta conservations, without significant increase in the 
calculation cost (see TableQ]in §I3.21 i. Note that too large val- 
ues of /, say / > 4, lead to the violation of the Courant con- 
dition. It is essential that the time-step of a particle j shrinks 



3 Recently, Springel 12009) has developed a mesh-based scheme employ- 
ing individual time-steps where the smaller time-step is adopted for integra- 
tions between neighboring meshes. This scheme does satisfy Newton's third 
law. 
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when the time-step of its neighbor particle i suddenly shrinks 
by a large factor. Thus particle i should let particle j respond 
to the change of its time-step. The key point of the limiter is 
to reduce effects of the violation of Newton's third law. Al- 
though our approach does not guarantee the third law, we can 
achieve sufficient accuracy by controlling the violation via the 
control parameter, /. 

We implement this limiter in the following way: when in- 
tegrated, particle i sends its new c/f, to its neighbor (y'-th) 
particles. The particles receiving the time-step of particle 
; compare it with their time-steps, dtj. When dtj > f dti 
they shorten their time-steps so that they satisfy the relation 
dtj =f dti. 

Note that this reduction of the time-step of particle j is only 
possible if the times of two particles and tj satisfy the con- 
dition 

ti>tj + fdt;. (3) 

If the above condition is not satisfied, the reduction of time- 
step results in the new time of particle j before the current 
time of particle i, requiring the backward integration of the 
entire system. 

In this case, the new time of particle j is set to a value that is 
consistent with the current system time (f,), and the time-step 
is set to the difference between particle' current time and new 
time. In Figure Q] we show schematic pictures of the tradi- 
tional individual time-steps method and our implementation 
of the individual time-steps method. 

3. NUMERICAL EXPERIMENTS 
3.1. Initial Setups and Method 

In order to verify the effect of the limiter, we performed 
two explosion tests. The first test deals with hydrodynamics 
only while the second one treats both hydrodynamics and self- 
gravity. For each test, we performed three runs: (a) integrated 
with global time-steps (hereafter Run a), (b) traditional indi- 
vidual time-steps (case 1. in Fig. [T] Run b), and (c) individual 
time-steps with the limiter (case 2. in Fig. [T] Run c). 

For the first test we adopted the Sedov problem, which is 
a pure hydrodynamical evolution of the point- like explosio n 
in the cold and homogeneous ambient medium dSedovll959l) . 
We prepared a glass distribution with 64 3 SPH particles. We 
enhanced the thermal energy of the central 32 particles in an 
SPH manner. The total thermal energy of the explosion is set 
to be unity. The internal energy of the ambient gas particle is 
set to be 10~ 6 of that for the hottest gas particle, although the 
original Sedov problem adopted a zero-energy background. 
Hence, the energy ratio is comparable with SN explosion, 
which generates a small amount of very hot gas (T ~ 10 7 K), 
in cold (T ~ 10 K) ambient gas. In this configuration, the po- 
sition of the shock front at the time T is ~ 1.15 x T 2 / 5 , when 
we assume the adiabatic index 7 = 5/3. 

In figure [2] we show the initial distribution of time-steps 
for the Sedov problem with individual time-steps as a func- 
tion of the radius. The difference in time-steps between the 
hottest particle and the ambient particles is ~ 10 3 . Note that 
time-steps of particles in the contact region surrounding the 
hot region are smaller than those of distant ambient particles 
(dt = 0.009), since we estimate these time-steps after the in- 
ternal energy of particles are specified. This distribution of 
time-steps is quite moderate compared to that in practical sim- 
ulations just after SN explosions, in which only the particles 
which receive the energy have small time-steps. 



We then perform an explosion test of a cold cloud including 
the self-gravity. For this test, we used the pa rticle distribut ion 
of the three-dimensional collapse test (e.g., E vrard|[T988l) at 
T = 3 where the system becomes the state of virial equilib- 
rium with the total energy of E ~ -0.6. We added the thermal 
energy to the central 32 particles in an SPH manner in order 
to drive the explosion of the system. The total energy of the 
system was set to be E = 10. 

We designed this test to mimic the SN energy input in high- 
density star-forming molecular cloud. The total binding en- 
ergy of a cloud with the t otal mass of 10 4 Mp, a nd the size 
of 10 pc (see Table 1 in Bergin & Tafallal 12007) is around 
-10 48 ergs, while the energy input of a type II SN is 10 51 ergs. 
Therefore, if we scale the energy accordingly, the energy in- 
put should be ~ 600. We used a much smaller value, so that 
the breakdown of the traditional method is less severe. We 
reset T = at the initial state. The number of SPH particles 
used here is 30976. The gravitational softening, e, is set to be 
0.05. 

We use the standard SPH method (e.g., lMonaghanlll992l) . 
The asymmetric form is adopted for the derivation of the in- 
ternal ^nergy^_The_signal velocity based artificial viscosity 
term (Monaghan 1997) is used. The viscous coefficient, a, is 
set to be two for the first test and unity for the second test. We 
determine the interaction (kernel) size of each SPH particle to 
keep the neighbor number t o 32 ±2. . Gravit y is solved by the 
Tree with GRAPE method (lMakinolll991bl) with an opening 
angle of 0.5. 

The time-integrator uses the leap-frog method. In runs of 
pure hydrodynamics, the time-step for individual particles is 
calculated as the smaller of the ones calculated using the sig- 
nal velocity based method (Mona ghanll 19971 : ISpringelll2005l) 
and th e acceleration time-step; 0.3 x (2/!/|flh|) ' 5 (Monaghan 
fl992h . Here h and flj, are the kernel size and the accelera- 
tion caused by the pressure gradient, respectively. The coef- 
ficient for the Courant condition, which is estimated by the 
signal velocity based method, is set to be 0.3. When we 
consider self-gravity, we modify the acceleration time-step 
taking into account the gravitational acceleration, a g , and e; 
0.3x(min(2/i,e)/|a h + a g |) - 5 . 

3.2. Results 

Figure [3] shows the result of the Sedov problem. In the 
top row, particle distributions of sliced regions (|z| < 0.01) 
are shown, and in the bottom row radial density profiles are 
shown. Run (a) (the left panels in Fig. [3]) shows a clear spher- 
ical structure and its density profile traces the analytic solution 
very well. 

Run (b) (the middle panels), however, displays penetrations 
of the hot gas particles into the ambient medium and fails to 
reproduce the analytic solution. This is because the time-steps 
of the gas particles composing the ambient medium are very 
long {dt ~ 10~ 2 ) compared with those of the initially heated 
particles (the minimum value is dt ~ 10~ 5 ) and therefore the 
ambient medium has had only two time-steps before the time 
of the snapshot (T = 0.02). They could not respond correctly 
to the expansion of the central hot region. 

Run (c) (the right panels) does not show this problem. In 
this run, the penetration of the hot gas particle occurring with 
the traditional method is suppressed perfectly and the ana- 
lytic solution of the Sedov problem is reproduced very well. 
This is because the short time-steps propagate to surround- 
ing particles quickly enough. As a result, cold gas particles 
can respond to the pressure of the hot gas particles. Conse- 
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FIG. 1. — Schematic pictures of individual time-step methods without and with the time-step limiter (case 1. and 2.). Particle (' suddenly shrinks its time-step 
(dtj) at the time indicated by the star symbol. In case 1., particle j does not respond to the change of dtj. In case 2., particles i and j always send their time-steps 
(dtj and dtj) to their neighbors when they become active. If the minimum time-step of neighbors is shorter than its own time-step by a factor / (see text), a 
particle shortens its time-step. 
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FIG. 2. — The distribution of initial time-steps for the Sedov problem with 
individual time-steps as a function of the radius. 



quently this run can reproduce the analytic solution of the Se- 
dov problem, like Run (a). We can conclude that this limiter 
or a similar criterion is necessary for the simulations of SPH 
with individual time-steps involving strong explosions. 

The total (kinetic + thermal) energy and linear momenta of 
the system at the final phase (T = 0.04) are shown in Table 
Q] Run (a) conserves the total energy to ~ 10~ 3 and the linear 
momentum with a very high precision. While Run (b) shows 
violations of conservation of both the total energy and the lin- 
ear momentum because the time-steps of cold gas particles 
are far too long to resolve the propagation of the blast wave. 

Introduction of the time-step limiter greatly reduces these 
errors. Runs (c) conserve the total energy to a level similar 
to that of Run (a). In addition, the conservations of the linear 
momentum in these runs are improved by four orders of mag- 
nitudes compared with that in run (b). When we compare runs 
with / = 2 and 4, we found the run with f =2 shows results 
slightly better in conservations. The distribution of particles 
and radial profiles are indistinguishable. Therefore we con- 
clude that the time-step limiter with / = 4 is acceptable. 

The calculation cost directly relates with the value of /; 
smaller / leads to larger calculation cost. Indeed, Run (c) with 
/ = 2 takes 1036 sec for integration until T = 0.4 whereas Run 
(c) with / = 4 requires only 489 sec. Interestingly, calculation 
costs for Runs (c) are smaller than that for Run (b). This is 
because the hot bubble does not expand correctly unless we 
introduce the limiter. 



TABLE 1 

Energy and momentum conservation and calculation cost. 



Run 




|(£i-£ F )/£i| 


Ipf|x 


\Pf\y 


\PfU 


Time (sec) 


a 




8.2e-4 


1.4e-16 


1.4e-15 


7.7e-16 


5646 


b 




7.2 


5.4e-l 


2.6e-l 


3.1e-l 


2068 


c(/ = 


2) 


7.6e-4 


3.2e-5 


3.1e-5 


l.le-5 


1036 


c(/ = 


4) 


5.9e-3 


1.7e-5 


2.7e-5 


2.2e-6 


489 



Note. Ei and Zip are the total energy at T = (initial phase) and 0.04 (final 
phase), respectively. |pf|x, IpfIy. and |pf|z indicate the absolute values of 
linear momenta for x, y, and z directions, respectively. 

Figure |4] shows sliced (|z| < 0.2) particle distributions for 
Runs (a), (b), and (c) (top to bottom) at selected epochs, for 
the second model with self-gravity. Since the total (kinetic + 
thermal + potential) energy of the system is positive and suf- 
ficiently large, the strong shock propagates outward against 
the self-gravity resulting in an expanding shell. Run (a) (top 
panels) shows this behaviors clearly. The behavior of the ex- 
pansion for Run (b) (middle panels) surprisingly differs from 
that for run (a). Again, the penetrations of the initial hot gas 
particles into the cold region are observed, and defects of the 
shell are shown in the middle-center and right panels. The 
introduction of the time-step limiter drastically improves the 
results as can be seen in the bottom panels of Figure|4] There 
are no penetrations or defects of the spherical shell-like shock 
which are visible in middle panels. 

These results clearly indicate that the use of individual 
time-steps without a time-step limiter is very dangerous for 
SPH simulations involving strong explosion problems, such 
as high-resolution SPH simulations of galaxy formation with 
a multiphase interstellar medium under 10 4 K. On the other 
hand, for traditional simulations of galaxy formation, the 
problem might not be too severe since these simulations adopt 
the cut off temperature of 10 4 K for cooling function, and tem- 
perature difference is at most around 100. We strongly rec- 
ommend that users of SPH who investigate systems in which 
"strong" explosions take place with individual time-steps im- 
plement some form of time-step limiter similar to what we 
introduced here. 
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FIG. 3. — Left, middle, and right column of panels show the evolution of particles integrated with (a) global time-steps, (b) individual time-steps, and (c) 
individual time-steps and the time-step limiter, respectively. The upper row shows particle distributions for three runs at T = 0.02. Points indicate projected 
positions of SPH particles in sliced regions (\z\ < 0.01) while crosses show those of initially heated SPH particles throughout whole regions. The lower row 
shows radial density profiles for three runs at T = 0.02. Points show geometrical means of densities for every radial bin with the width of 0.005 from the origin 
of coordinates. Solid (gray) curves are the analytic solution for the zero-energy background obtained by Se dovl (19591) . 



putations were carried out on XT-4/GRAPE systems (project CfCA, of the National Astronomical Observatory of Japan. 
ID:g08al9) at the Center for Computational Astrophysics, 
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FIG. 4. — Evolution of particles for three runs. Upper, middle, and lower 
rows of panels show the evolution of particles integrated with (a) global time- 
steps, (b) individual time-steps, and (c) individual time-steps and the time- 
step limiter, respectively. Points and crosses indicate projected positions of 
SPH particles in sliced regions (|z| < 0.2) and those of initially heated SPH 
particles in whole regions, respectively. 



